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ABSTRACT 

We present an analysis of the level of polarized dust and synchrotron emission using the WMAP9 and 
Planck data. The primary goal of this study is to inform the assessment of foreground contamination 
in the cosmic microwave background (CMB) measurements below £ ^ 200 from 23 to 353 GHz. We 
compute angular power spectra as a function of sky cut based on the Planck 353 GHz polarization 
maps. Our primary findings are the following. (1) There is a spatial correlation between the dust 
emission as measured by Planck at 353 GHz and the synchrotron emission as measured by WMAP 
at 23 GHz with p ^ 0.4 or greater for ^ < 20 and /sky > 0.5, dropping to p « 0.2 for 30 < ^ < 200. 
(2) A simple foreground model with dust, synchrotron, and their correlation fits well to all possible 
cross spectra formed with the WMAP and Planck 353 GHz data given the current uncertainties. (3) 
In the 50% cleanest region of the polarized dust map, the ratio of synchrotron to dust amplitudes at 
90 GHz for 50 < <110 is 0. 3 / 0 . 2 • Smaller regions of sky can be cleaner although the uncertainties 

in our knowledge of synchrotron emission are larger. A high-sensitivity measurement of synchrotron 
below 90 GHz will be important for understanding all the components of foreground emission near 90 
GHz. 


1. INTRODUCTION 

The drive to characterize the GMB temperature and 
polarization anisotropy has led to a steady increase in 
our knowledge of foreground emission from our galaxy. 
In this paper we focus on the diffuse polarized fore¬ 
ground emission in the 23-353 GHz regime as measured 
by the WMAP (Bennett et al.2013, hereafter B13) and 
the Planck (Planck Gollaboration et al. 2015a) satel¬ 
lites. There are two well established polarized compo¬ 
nents: synchrotron emission and thermal dust emission. 
It is known that they are spatially correlated as shown 
in Kogut et al. (2007), Page et al. (2007), and Planck 
Gollaboration et al. (2015c) (hereafter K07, P07, and 
P15), that the polarized synchrotron is described by a 
power law with the spectral index varying over differ¬ 
ent galactic latitudes (K07; B13), and that the polarized 
thermal dust emission is well-described by a blackbody 
emitter with an emissivity that scales as a power law in 
frequency. 
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with Tdust = 19-b K and pd = 1.59 (P15), where Bj^{T) is 
the Planck function and = 353 GHz is the reference 
frequency. 

We still have a lot to learn about foreground emis¬ 
sion. For example, there is no definitive explanation for 
13d = 1.59 or for the observation that the spectral index 
of polarized synchrotron is different than that of unpo¬ 
larized synchrotron (Gold et al. 2011). It is also seen that 
the power in EE is roughly twice that in BB in both syn¬ 
chrotron and dust (P07; Gold et al. 2011; Planck Gollab¬ 
oration et al. 2015b; Planck Gollaboration et al. 2014b, 
hereafter P14). These observations are likely rooted in 
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the configuration of the galactic magnetic field relative 
to the galactic plane and the coherence length of the 
magnetic field relative to the dust and electron distribu¬ 
tions. Recently, for example, Planck Gollaboration et al. 
(2015d) showed that the EE/BB ratio for dust may arise 
from the correlation between the filaments and the mag¬ 
netic field orientations. Nevertheless, our characteriza¬ 
tion of foreground emission is largely empirical. 

One of the notable findings is shown in Eigure 16 of 
B13. The “synchrotron spectrum” appears to flatten as 
one goes to higher frequencies and there is an apparent 
dust contribution. With synchrotron and dust parame¬ 
terized as Tg (X and Td oc a possible model is 
put forward that combines the Ps flattening in frequency 
from —3.13 at 33 GHz to —2.43 at 94 GHz with an ad¬ 
ditional dust component that rises with = 1.4 — 1.5. 
The WMAP team cautiously interpreted this flattening 
as the synchrotron template tracing some of the dust due 
to shortcomings of the dust template. With the Planck 
353 GHz polarization maps, we can now better test this 
model. 

Our main goal is to test a simple model of polar¬ 
ized dust and synchrotron emission over different frac¬ 
tions of sky and angular scales. Often foreground stud¬ 
ies are done with templates in map space. This proce¬ 
dure weights the largest angular scales because they have 
the highest signal to noise ratio. Here we examine the 
emission characteristics in £ ranges for masks of the sky 
relevant to searches for primordial B-mode polarization. 

2. DATA 

WMAP delivered all-sky polarization maps observed 
with ten independent differencing assemblies (DAs) cov¬ 
ering five different frequency bands. We use the single 
year DA maps to examine the map statistics. We And 
this important for a thorough assessment of the angu¬ 
lar power spectrum uncertainty. In particular, it gives 
us insight into the cross term between the signal and 
noise, and also allows us to investigate the distribution 
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Fig. 1.— Masks showing regions with incrementally larger levels 
of polarized emission in the Planck 353 GHz maps, with /g^y of 0.3, 
0.5, and 0.7 (see details in Section 2). These regions do not neces¬ 
sarily indicate progressively higher levels of polarized emission for 
50 < ^ < 110. We also show the approximate regions observed by 
the BICEP2 and Keck experiments (Ade et al. 2015) and SPIDER 
(Fraisse et al. 2013). 

of the noise. Planck has released all sky polarization 
maps at 30, 44, and 70 GHz from the Low Frequency 
Instrument (LFI) and 100, 143, 217, and 353 GHz from 
the High Frequency Instrument (HFI). For this work, we 
use the full-mission and half-mission splits for only the 
353 GHz maps. Intensity to polarization leakage in the 
353 GHz maps, if not corrected, corresponds to < 4% 
of the EE/BB power spectra for > 50 in the similar 
sky regions to what we define below (P14). We use the 
corrected maps. The precise level of the residual leakage 
at large angular scales is not yet published by the Planck 
collaboration. 

We make masks with the Planck 353 GHz polarization 
map^, P = y/Q‘^ P obtained from first degrading 

the Stokes Q and U maps to HEALPix (Gorski et al. 
2005) A^side = 128 then smoothing them with a Gaus¬ 
sian EWHM of 1°. We retain pixels above different 
thresholds set to get masks^ with /sky of 0.7, 0.5, and 
0.3 as shown in Eigure 1. The regions defined this way 
sequentially mask more of the galactic plane and high 
polarized emission regions. This approach selects the 
regions with incremental levels of power in the large 
angular scales, as these modes have the highest signal 
to noise ratio. We apodize the three masks with a 
Gaussian EWHM of 5°. Since we are interested in 
analyzing the diffuse foreground emission, we mask all 
sources found in the WMAP point source catalogue 
and sources with S/N >3 and flux density above 
400 mJy from Planck catalogues of Gompact Sources 
and Galactic Gold Glumps at 353 GHz within a radius 
of 15', then apodize the edges smoothly to unity with a 
Gaussian a of 0.559°. This particular a ensures that the 
argument of the exponential is at —5 by a radius of 1.5° 
from each source. We normalize the resulting mask then 
compute /g]^ from /skyii’ 2/'^4 (equation (17) in Hivon 
et al. 2002), where Wi is the ith moment of the mask, 

^ We do not subtract the noise bias and the CMB from P and 
this may have an effect on the exact definition of the regions. 

^ We first make a binary mask by padding pixels over some 
threshold with —1 and the rest with 1, then smooth the mask, 
set the negative pixels to 0 and retain the positive pixels to form 
largely contiguous masks with /sky of 0.7, 0.5, and 0.3. 


to find 0.62, 0.45, and 0.28 for our three masks. The 
main results presented in the following sections are not 
sensitive to the exact details of the apodization used here. 

3. METHOD 

We use the nominal full sky pseudo-Q approach 
(Hivon et al. 2002; Kogut et al. 2003; Brown et al. 2005), 
which accounts for incomplete sky coverage and beam 
smoothing to give unbiased estimates of the angular 
power spectrum. We confirmed the accuracy of our code^ 
with an independent code from Kendrick Smith, and aiso 
reproduced the 353 GHz power spectra shown in Eigure 
2 in P14 up to smaii differences seemingiy due to the 
differences in the assumed masks. The pure pseudo-C^ 
estimator is not needed in this anaiysis since the power 
in E- and B-modes are simiiar for foregrounds and our 
masks retain iarge sky fractions. 

We use only cross spectra in estimating the signal to 
avoid noise bias in our estimates. The power spectrum at 
353 GHz is obtained from the cross spectra of the half¬ 
mission maps. Eor the cross spectra between WMAP 
and Planck, we use the WMAP single year maps and the 
Planck half-mission maps, then give the weighted mean 
for each cross-frequency spectra, accounting for the co- 
variance among many cross spectra as we describe below. 
Eor all auto-^ and cross-DA spectra within WMAP, we 
use the nine single year maps to get 36 cross spectra for 
each DA, and 81 cross spectra for each cross-DA pair. 
Then we report the five single- and ten cross-frequency 
spectra for WMAP as the weighted mean of all spectra 
involving the relevant DAs. The W4 channel is omitted 
due to its failure of internal consistency tests (Jarosik 
et al. 2011), which we independently confirmed. This 
channel had a large in-flight 1 //knee compared to those of 
the other DAs (Jarosik et al. 2003). We also note that all 
maps are used at their provided resolution: A^side = 512 
for WMAP and A^side = 2048 for Planck. 

Our goal is to parametrize the foreground emission, a 
non-Gaussian signal, with representative error bars. To 
this end we make a number of approximations in esti¬ 
mating the statistical uncertainty as described below. In 
particular, for WMAP we derive the statistical uncer¬ 
tainty from the distribution of the spectra formed with 
the nine single year maps, as each measures essentially 
the same sky (in one year, WMAP measures the full sky 
twice). A similar treatment is shown in Hinshaw et al. 
(2009), although it is for / < 10 and for foreground- 
cleaned maps, which should be close to Gaussian. De¬ 
tailed studies at large angular scales, say < 20, need 
to use the pixel-based likelihood solution with the pixel- 
based noise covariance matrix (Appendix D of P07) or 
even better, simulations based on the time-ordered data. 
However, for our purposes, the pixel-based likelihood is 
not needed. 

To increase the signal to noise ratio, we compute the 
power spectra binned in multipole bands b of width 
specified below. The auto power spectrum of a single 

The code uses HEALPix libraries for spherical harmonic trans¬ 
forms and SLATEC Fortran subroutine, DRC3JJ.F9, for the com¬ 
putation of Wigner-3j symbols. 

^ We use the term auto-DA spectra to refer to the cross spectra 
associated with one DA. The term auto power spectrum refers to 
the power spectrum of a map of signal plus noise. 
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Fig. 2.— The beam-deconvolved noise power spectra of the 9- 
year K-band maps are shown for the mask with of 0.62. The 

solid red line corresponds to from subtracting each C}) 

estimated from the 36 cross spectra of the single year maps from 
each auto power spectrum of the 9-year maps. The black dashed 
line corresponds to obtained from the mean of each 

spectra from 1000 MC simulations of the noise maps. A similar 
quantity illustrating the large scale noise correlation is shown in 
Figure 16 of POT. 

map in bin b gives (7^ + where is the binned signal 
power spectrum, which includes both foregrounds and 
the CMB, and is the binned noise power spectrum. 
We can estimate Ni, from the auto power spectrum by 
subtracting Cb estimated from the cross spectra. The 
estimate of here is then sensitive to any 1/^-like com¬ 
ponent of the noise, since the particular realization of 
the noise contained in the data retains the large scale 
correlations from the 1/f noise in the time-ordered data. 
Alternatively, we can estimate N^, from the Monte Carlo 
(MC) simulations of the noise sky using the provided 
noise maps (the diagonal element of the full noise covari¬ 
ance matrix). These simulations, without the full noise 
covariance matrix, do not capture the l/£ component of 
the noise, as shown in Figure 2 for the 9-year K-band 
maps. 

To determine the uncertainty for a set of cross spectra 
from many data splits, we need the covariance matrix 
that describes the correlations among the cross spectra 
with common data splits. The usual Gaussian approxi¬ 
mation of the total uncertainty for an auto power spec¬ 
trum in bin b is given by 

(ACbf = -iCb+Ni,f, ( 2 ) 

with Vb defined as. 


Ub = (24 + 1)^4/s®]^5 

where £b is the central multipole in 6, and is the 
effective fraction of sky as described in Section 2. The 
term is the sample variance. We generalize this ex¬ 
pression for cross spectra uncertainties following the ap¬ 
pendix in Das et al. (2011). The covariance matrix ele¬ 
ment of cross spectra for each bin b is given by 

^ijBxlD)\ 

(3) 


^{iAxjB)-{kCxlD) _ 
^bb ~ 

Uh 


y^(^iA X kC^ 


where different instrument channels and data splits are 
labeled with uppercase and lowercase roman letters re¬ 
spectively. Note that we recover equation (2) for the case 
of auto spectra, i = j = k = I and A = B = C = D, 
since ^(^bAxzA)^ = Cb ^ Nb. 

When evaluating the distribution of the cross spectra 
from multiple splits of data, there is no sample variance 
as each data split contains the same signal spectrum; 
however, there are cross terms CbNb between signal and 
noise, different in each spectrum. We illustrate this by 
considering equation (3) for WMAP K-band for bin 6, a 
36 x36 covariance matrix (for the spectra with nine single 
year maps excluding the auto spectra). The variance of 
year 1x2 is 


^{2KxlK)-{2KxlK) 


1 

T^b 


^(2Kx2K)\^^^(lKxlK) 


(4) 


which evaluates to 


^{2KxlK)-{2KxlK) 

^bb 




+ Nl^) 


Similarly, the covariance of year 1x2 with year 1 x 3 is 


^{2KxlK)-{3KxlK) 

^bb 



We can see that the matrix has the form where the diag¬ 
onal elements include the N x N and the C x N terms, 
and the off diagonal elements have just the C x N terms 
or zeros. 

The effective number of modes z/^, parametrized by 
our estimate of needs to be calibrated with sim¬ 

ulations to account for the spatial inhomogeneity and 
correlations in the noise and the complex mask. In par¬ 
ticular, to calibrate Vb a-f all angular scales one needs the 
pixel-based covariance matrix to simulate the noise maps 
and a template for the signal. With these, the distribu¬ 
tion of the power spectra has the representative level of 
variance from the C x N and the N x N terms, but 
does not have sample variance. For the native resolution 
and our £ range, such a calculation is computationally 
intractable. Hence we instead use the uncertainties asso¬ 
ciated with each pixel (the diagonal elements of the full 
pixel-based covariance matrix) including the QU corre¬ 
lations to generate the noise map and simulate the signal 
sky from a power law in £ fitted to the estimates of 
Cb for each single- and cross-frequency spectra. Then we 
get the statistical variance of the power spectrum of the 
simulations, (AC^^)^, by subtracting the sample vari¬ 
ance (obtained from the signal-only simulations) from 
the variance in the signal plus noise simulations. The 
/g®^ calibration factor is given by the ratio of the statisti¬ 
cal variance in the simulations to the analytic expression 
from equation (3), 

9b — [^^b J /'^hh ’ 

where the signal and noise power spectra in 

and are from 2000 simula¬ 

tions. Note this method of obtaining gb takes into 
account the effects from the noise inhomogeneity and 
the QU correlations but not the l/£ noise component. 
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TABLE 1 

Single and cross-frequency spectra x^/dof 


Channel 

u [GHz] a 

dof*^ 

i e [6, 23] 

£ e [24,49] 

i e [50, 110 ] 




EE 

BB 

EE 

BB 

EE 

BB 

K 

22.4 

35 

2 . 10 ^ 

2 . 75 ^ 

1.64 


1.77 


1.43 


1.20 


0.77 


0.58 


0.99 


0.78 


1.11 


0.91 

KKa 

27.1 

80 

6.24 

7.29 

4.79 


5.09 


1.52 


1.13 


1.39 


0.98 


1.09 


0.74 


2.02 


1.43 

KQ 

30.1 

161 

8.15 

9.02 

2.23 


2.40 


2.06 


1.51 


1.97 


1.34 


1.14 


0.73 


1.27 


0.83 

KV 

36.7 

161 

4.96 

5.33 

2.64 


2.84 


1.56 


1.24 


1.75 


1.23 


1.11 


0.68 


1.33 


0.84 

KW 

45.8 

242 

4.15 

4.34 

3.04 


3.30 


2.18 


1.57 


1.42 


0.94 


1.31 


0.83 


1.33 


0.81 

Ka 

32.7 

35 

1.31 

1.40 

1.91 


1.60 


0.72 


0.55 


1.32 


0.94 


0.80 


0.65 


1.20 


1.00 

KaQ 

36.3 

161 

2.07 

1.85 

1.30 


0.99 


1.40 


1.02 


1.28 


0.88 


1.45 


1.04 


1.28 


0.99 

KaV 

44.3 

161 

1.63 

1.40 

1.41 


1.10 


1.35 


0.92 


1.36 


0.97 


1.26 


0.91 


1.26 


0.91 

KaW 

55.2 

242 

1.69 

1.45 

1.53 


1.14 


1.42 


1.06 


1.45 


1.02 


1.21 


0.83 


1.37 


0.95 

Q 

40.3 

152 

2.01 

1.40 

1.82 


1.24 


1.39 


0.96 


1.56 


1.08 


1.47 


1.13 


1.29 


1.01 

QV 

49.3 

323 

1.89 

1.27 

1.64 


1.11 


1.40 


0.96 


1.52 


1.01 


1.09 


0.83 


1.42 


1.08 

QW 

61.4 

485 

1.74 

1.16 

2.02 


1.27 


1.51 


1.04 


1.36 


0.89 


1.38 


1.04 


1.26 


0.94 

V 

60.2 

152 

1.61 

1.05 

2.20 


1.43 


1.50 


1.02 


1.35 


0.91 


1.23 


0.99 


1.39 


1.09 

vw 

75.0 

485 

1.77 

1.14 

2.43 


1.56 


1.42 


1.02 


1.73 


1.27 


1.21 


0.92 


1.21 


0.91 

w 

93.4 

350 

1.63 

1.05 

2.95 


1.93 


1.48 


1.01 


1.52 


1.08 


1.28 


0.99 


1.09 


0.85 

K353 

90.4 

17 

1.10 

1.14 

1.01 


1.32 


1.22 


1.06 


1.39 


1.25 


0.87 


0.77 


0.90 


0.73 

Ka353 

109.1 

17 

1.51 

1.66 

1.69 


1.88 


0.44 


0.38 


0.65 


0.55 


1.15 


1.08 


0.51 


0.46 

Q353 

121.2 

35 

1.08 

1.21 

1.15 


1.28 


1.49 


1.34 


0.96 


0.86 


1.02 


0.98 


0.61 


0.55 

V353 

148.1 

35 

1.36 

1.52 

1.08 


1.15 


0.88 


0.81 


1.40 


1.27 


0.64 


0.59 


0.85 


0.75 

W353 

184.4 

53 

1.67 

1.79 

0.98 


1.08 


1.27 


1.15 


1.52 


1.41 


1.09 


1.06 


0.91 


0.84 

353 

364.2 

0 























^The effective frequency centers for the given foreground spectrum (see details in Section 3). The single-frequency spectra show the 
obtained physical frequency centers and each cross-frequency spectrum shows the geometric mean of the two frequencies for the pair. 
^The degrees of freedom (dof) for the x^ in each set of spectra, accounting for the removal of the mean. 

^The reduced x^ for single and cross-frequency spectra in the of 0.62 region. Boldface denotes > 10% difference between x^/dof 
from using the full covariance matrix and using just the diagonal elements. 

Shaded columns are x^/dof after calibration from gi). See Section 3 for details. 


For one ^-bin 6 , there are 20 covariance matrices^ of di¬ 
mension dof+ 1 , where dof is given in the third column of 
Table 1. We construct each covariance with the estimate 
of initially from the mean of all spectra involving 
the relevant DAs, and the estimate of Ni, from the auto 
spectra, as described above (thus here contains the 
l/£ component of the noise). Then with each auto-DA 
covariance matrix, we get the weighted mean of the cross 
spectra to estimate use this estimate to rebuild the 
covariance matrix, then iterate until convergence^. This 
second step typically shifts C 5 by less than O.lcr, though 
the shift could be as much as 2.8cr. The shift essen¬ 
tially shows the difference between the weighted mean 
and the initial arithmetic mean of the cross spectra. The 
final estimates of the auto-DA C 5 obtained this way are 
the signal input for the cross-frequency covariance ma¬ 
trices. When the WMAP channels are noise-dominated, 
the auto-DA Cb can fluctuate negative. In such case, we 
simply use the auto spectra as the estimate of and also 
ignore the Cb x Nb terms when building the covariance 
matrices to avoid any non-physical noise representation 
in the data. This happens in 53 out of 162 auto-DA Cb 
for all bins and masks^, with 31 of them in the W-band 
and the rest in Q or V bands. 

The bandpass mismatch within each band, say, be¬ 
tween the WMAP channels Q1 and Q 2 , leads to a small 
amount of fluctuation among the single-frequency spec¬ 
tra and the cross-frequency spectra with the other chan¬ 
nels. WMAP calibrates on the CMB. However, the low 
frequency channels are dominated by synchrotron emis¬ 
sion and have different central frequencies. Thus, Q 1 
and Q 2 have slightly different synchrotron levels. We 

® For all single- and cross-frequency spectra except 353 GHz. 

^ Convergence is achieved with 3 iterations. 

8 9 (DA) X 2 (EE, BB) x 3 (Gbins) x 3 (masks) = 162 


have combined all of the cross spectra from Ql, Q 2 , 
and QlxQ 2 , then give the effective central frequency in 
Table 1. We do not scale the spectra from the differ¬ 
ent channels before combining them, and this affects our 
X^/dof by < 1%. The frequency centers in Table 1 were 
determined iteratively considering the spectral index of 
the signal in antenna temperature^. For WMAP, we use 
the software on LAMBDA^^ to get the corrections, which 
are based on Jarosik et al. (2003). For Planck, we use 
the corrections in Table 4 of Planck Collaboration et al. 
(2014a). This ensures our estimates are the representa¬ 
tive power at the given effective central frequencies. 

We use the covariance matrices with no sample vari¬ 
ance to evaluate the consistency of the measurements and 
the error bars. We show the reduced for all single- and 
cross-frequency spectra in Table 1 for the mask with 
of 0.62 before and after calibrating z /5 with gb. The cross¬ 
frequency spectra here combine all corresponding cross- 
DA spectra. The left hand column in each pair shows 
from before calibrating Ub. After correcting Ub for the 
number of modes that are lost, the scatter in the data at 
low i is large even though the C x N terms and the l/l 
component of the noise are accounted for. This is likely 
due to not using the full pixel-based covariance matrix. 
The departure from unity in the y^/dof corresponds to 
an underestimate of the uncertainty up to a factor of 3 
for WMAP and an overestimate up to a factor of 1.5 
in WMAP X Planck. In general, though, we find that we 
can treat the signal as normally distributed for i > 23. 

^ We use all spectra in the first Gbin for EE and BB simul¬ 
taneously to find the spectral index at each band, compute the 
frequency center corrections for the given indices, recompute the 
spectral indices at the new frequency centers, and iterate the pro¬ 
cess until convergence. 

http://lambda.gsfc.nasa.gov/product/map/ 
current/effective_freq.cfm 
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As a consistency check of the uncertainties, we com¬ 
pare our estimates from the data split maps (single year 
and half-mission maps) to those estimated with the full 
coadded maps (9-year and full-mission maps). First we 
subtract C 5 estimated from the data split maps from the 
auto spectra of the full coadded maps to get N}j. We then 
get the statistical uncertainty with equation (3) . For the 
cross-frequency spectra, the agreement of the uncertain¬ 
ties is typically within a few percent with a deviation of 
^ 30% occurring at Kx353 for the first i-hin in EE for 
each mask. Eor the single-frequency spectra, our nom¬ 
inal uncertainties using the data split maps are greater 
than those estimated from the full coadded maps by typ¬ 
ically ^ 5%, as the auto spectra are not included for the 
single-frequency spectra, while from the full coadded 
maps include the full data^^. 

Eor fitting to the model, described in the next section, 
we need the covariance matrix of the 21 spectra in each £- 
bin b. To form this we use the Nij from the full coadded 
maps, apply the ^5 factors, then scale the components 
appropriately to account for not including the auto power 
spectrum in our estimates of 

4. MODEL 

We seek a simple model to characterize the polarized 
foreground emission and to serve as a guide to the degree 
of cleaning necessary to identify primordial or lensing B- 
modes at < 150. A number of studies have been done 
previously for modeling and cleaning of the foreground 
(e.g., B13, P14, Euskeland et al. 2014, Dunkley et al. 
2009, Planck Collaboration et al. 2015b, and Watts et al. 
2015), to which we add by basing our model on the mea¬ 
sured set of all single- and cross-frequency spectra with 
Planck 353 GHz and WMAP, as a function of Abin b. 
One advantage of this approach is that variations with 
angular scale (or spectral bins) may be identified. 

At the current level of precision, synchrotron and dust 
emission may be modeled with the WMAP K-band and 
the Planck 353 GHz channels as 

Mi,,) = S{iy) + Diiy) 

= as{jy/jyK)‘^‘S + ad£d{i^,/3d)^, 

where M{u) is a Stokes Q or U map in antenna tem¬ 
perature at frequency u in GHz, is the Planck 

dust model from equation ( 1 ), S and D are normalized 
synchrotron and dust spatial templates at frequencies uk 
and 353 GHz, with amplitudes ag and ad- When taking 
the cross spectra between maps M(z^i) and M{y 2 )^ we are 
effectively looking at the covariance between the spatial 
templates weighted by the amplitude and the frequency 
dependence terms, 

^ al{viV2lvKY‘ + a\£d{vi, Pd)£d{v2, Pd) 

+PsdCisCid{{vi/V k)^‘£ d{v2, Pd) + {^2/Vk)^‘£ d{vi, Pd)), 

( 6 ) 

where {ad) is defined as ag (ad) times the variance in 

For instance, out of the total 45 possible spectra within 
WMAP K-band, 9 auto spectra are not included. This amounts to 
a factor of a/ 9/8 larger error bar for the case of noise domination, 
compared to that from using the full 9-year map. 


§ (B), and pgd is defined as, 

Psd =< S X D > /v< S X S >< D X D >, ( 7 ) 

the covariance between dust and synchrotron templates 
divided by the square root of the product of the variance 
of each (free of noise bias). However, when we fit to this 
model, Pgd is a free parameter that we check with the 
definition. 

As shown in the next section, we can fix Pd = 1-59 
(P15) and obtain good fits between all single- and cross¬ 
frequency spectra formed with the WMAP and the 
Planck 353 GHz maps. Thus we consider just four pa¬ 
rameters (Sg, yds, ad^ psd) as our nominal model, which 
we denote model A. Alternatively, we can fix pgd from 
equation (7) using the WMAP K-band and the Planck 
353 GHz maps as the synchrotron and dust templates, 
and use (S^, yd^, Pd) as the free parameters, which 
we denote model B. We do not need more than four pa¬ 
rameters for the 21 single- and cross-frequency spectra 
in each Gbin. 

Before fitting to the model, we subtract lensed ACDM 
GMB power with the tensor to scalar ratio r = 0 from 
EE and BB, and add the appropriate cosmic variance 
to the uncertainty. We then convert these foreground- 
only spectra into antenna temperature units, and fit the 
resulting set of spectra for each mask to equation (6), 
using the full covariance matrix for each bin (described 
in the end of Section 3). 

The model parameters and the uncertainties are ob¬ 
tained from two ways. Eirst, we simulate the data 10000 
times using the covariance matrix, then for each simula¬ 
tion the parameters are found by minimizing the 
the fit. We then give the central values and uncertain¬ 
ties from the maximum and Icr of the distributions for 
each parameter. The priors we adopt for this method 
are uniform and largely unbounded except where they 
are non-physical, 

as>0 
> 0 
<0 
Pd>0 

\psd\ < 1 . 

We find consistent parameters from the posterior prob¬ 
ability distributions sampled with the Markov Ghain 
Monte Garlo (MGMG) methods using the EMCEE code 
(Eoreman-Mackey et al. 2013) with a conservative uni¬ 
form prior distribution of — 4 < < 0. 

5. RESULTS 

The power spectra in where Vi = i{£-\- l)Q/(27r), 
and the corresponding fits to model A (from mini¬ 
mization) are shown in Eigures 3, 4, and 5 for the three 
masks. The spectra are GMB-subtracted and shown in 
antenna temperature units. The panels under each fig¬ 
ure show the residuals on a linear scale. The square root 
of the diagonal elements of the covariance matrices are 
shown as the error bars. The single-frequency spectra are 
shown in red diamonds, and the cross-frequency spectra 
are in green circles, plotted at the geometric mean of 
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Fig. 3.— The data and model for = 0.62. All y-axes are 

power in in antenna temperature. The upper panel of each 

pair has a logarithmic y-axis whereas the lower panel is linear. 
Left and right columns correspond to EE and BB respectively. See 
Section 5 for details. 

the two frequencies. The cross spectra do not necessar¬ 
ily represent the power at these plotted frequencies. The 
error bars that hit the x-axes pass through the zero line, 
and the yellow error bars denote two sigma error bars for 
when the spectrum estimate is negative. 

In each panel, the blue line shows the ACDM CMB 
spectrum in antenna temperature. For BB, the magenta 
line shows the leasing B-mode signal plus an r = 0.1 
signal. The negatively-sloped black dashed line shows 
the synchrotron spectrum from 5^ and /3s . Similarly, the 
positively-sloped black dashed line shows the dust spec¬ 
trum from and pd = 1.59. The dust-synchrotron cor¬ 
related component is shown as a red dashed line. The ex¬ 
pected total foreground power at any frequency is given 
by the smooth solid black curve (the sum of the three 
dashed lines), from equation 6 for = z^ 2 . Hence the 
cross-frequency spectra (green circles) are not generally 
expected to lie on this line. The fitted model is shown 
by the zigzag cyan line with cyan open circles, capturing 
both the single- and cross-frequency spectra. We em¬ 
phasize the model plotted is valid only at the frequencies 
given in Table 1 indicated by the open circles, and should 
not be interpolated between them. Lastly, the residuals 
to the best-fit model (the zigzag cyan line with cyan open 
circles) are shown in the panels below the spectra. A sig¬ 
nificant level of constraining power is seen in the cross 
spectra between Planck 353 GHz and WMAP (five green 
circles between 90-180 GHz). 

The parameters of the best fit model found from 
minimization and MGMG are given in Table 2. The pa- 


Fig. 4.— The data and model for = 0.45 following the 

conventions of Figure 3. 
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Fig. 5.— The data and model for = 0.28 following the 

conventions of Figure 3. 
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TABLE 2 

Results for model A 
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^In antenna temperature units at the effective frequency center for K-band. 

^In antenna temperature units at the effective frequency center for 353 GHz. 

‘^Note the results and the probability to exceed (PTE) values for the different masks are not independent as they include common regions. 
We assume the same dof = 17 for the two methods. 

•^Unshaded columns show the parameters found from minimization. 

®Shaded columns show the parameters found from MCMC. 


rameters obtained the two ways are consistent. Table 2 
shows that this simple model is a good fit to the data. 
Figure 6 shows the fitted values of psd^ and Pd- The 
top panel suggests that /3s becomes more positive at finer 
angular scales although the large uncertainty precludes 
drawing a definitive conclusion. However, Ps = —3.0 
does not appear to be less likely than /3s = —3.3, espe¬ 
cially for 50 < ^ < 110. The center panel shows psd from 
model A and also obtained from the templates (equa¬ 
tion (7)) for EE and BB, showing consistency between 
the two methods. We find psd ~ 0.4 for 6 < ^ < 23 in 
the masks with of 0.62 and 0.45, then it drops to 
a constant value 0.2. While the particular region we 
consider for of 0.28 shows smaller psd in the first 
Abin than those in the other two masks, we find a re¬ 
gion with similar based on the dust intensity with 
Planck 857 GHz map also shows psd ~ 0.4. Note the 
largest angular scales not presented here show a larger 
level of correlation with psd > 0.5. The bottom panel 
shows the results for model B, with l3d freed and psd 
fixed from equation (7). We find consistent indices with 
those from P15. 

The 50 < < 110 region is the easiest place to search 

for primordial B-modes from the ground. We show the 
marginalized distributions on the parameters for of 
0.28 in this i range in Eigure 7. The current uncertainty 
on ^s foi* such is large for this i range, and this has 
a profound impact on the level of synchrotron emission 
near 100 GHz. A difference of 0.5cr in /3s can lead to a 
factor of ~ 4 times larger synchrotron power at 100 GHz. 
Eigure 8 shows the best fit spectrum to model A (from 
MGMC) for the three different sky cuts. Although the 
fits to the model were obtained with /3d fixed, we have 
also included the appropriate uncertainty on /3d for the 
given /gj^ (P14) in the shown uncertainty bands. The 


dark (light) gray bands correspond to la {2a) spread at 
each frequency in the individual fits from the accepted 
sets of parameters, illustrating the current uncertainty 
on polarized foregrounds in these regions. Euture stud¬ 
ies including the rest of the Planck polarization data will 
shrink these error bands. The red points in each panel 
show the best fit synchrotron and dust power at the refer¬ 
ence frequencies^^. The minimum frequency in the given 
i range for the three masks is ~ 75 GHz, similar to that 
for the temperature anisotropy (B13; Planck Collabora¬ 
tion et al. 2015b). Eor /g®j^ of 0.62 and 0.45 the ratio of 
synchrotron to dust emission at 90 GHz in amplitude is 
O.fi^o'i arid 0-3lo!2 respectively. In the of 0.28 mask, 
the uncertainty on the synchrotron emission is large. The 
ratio of synchrotron to dust amplitudes at 90 GHz is < 5 
at 95% confidence, and the ratio of the correlated com¬ 
ponent to the dust amplitudes is 0.21 q;}. 

The cyan lines in the top panels mark 90 and 150 GHz, 
typical frequencies for ground based observations. The 
index of the total foreground spectrum between these 
two frequencies, /dfg, is shown in the bottom panels for 
each mask. Blue vertical lines at Pfg ~ —0.7 indicate the 
index of the CMB spectrum between 90 and 150 GHz. 

Lastly, we examine the model in the region observed 
by the BICEP2 and Keck experiments (Ade et al. 2015) 
which is cleaner in dust in the 50 < ^ < 110 range than 
our /gj^ = 0.28 region. We define the region by set¬ 
ting the pixels inside the published field outline^^ to 1 
and the outside to —0.01 at A^side = 612, smoothing 
with 6° Gaussian EWHM, then setting negative pixels 

The frequency centers given for K-band and 353 GHz in Ta¬ 
ble 1. 

http://bicepkeck.org/B2_3yr_373sqdeg_field_ 
20140509.txt 



















Fig. 6. — The three panels show /3s, psd^ and fitted for the 
two models for the three masks in the three ^-bins. The results for 
each ^-bin are separated by the solid black vertical lines. For each i- 
bin, to the left (right) of each dotted black line are the parameters 
for EE (BB), where different colors represent the results for the 
different masks. In each groups of points, the squares (triangles) 
show the results from model A (model B) from minimization. 
We show the parameters for model A from MCMC for the lower 
signal to noise ratio regime in stars. In the of 0.28 region, /3s 
in the second ^-bin is a 95% confidence limit from the prior bound 
at —4. For model A, ^d is fixed, and thus not shown in the bottom 
panel. 


to zero. The of the mask is 0.0088. Before com¬ 
puting the power spectra in the region, all maps are 
high pass filtered in harmonic space with an exponen¬ 
tial function smoothly going from zero to one from I of 
2 to 21. The noise in WMAP is too large to make a 
significant measurement of in a single i-hin. There¬ 
fore we fix psd = 0.2, Pd = 1.59, and impose a Gaus¬ 
sian prior on Ps = —3.0T0.3 based on the results from 
the larger sky cuts, then fit the model using MCMC si¬ 
multaneously in 3 even i bins between 21 and 125 us¬ 
ing = A^{£/S0)~^-^^ (P14) to get the synchrotron 

and dust power at = 80 {As and Ad). The power 
law scalings is consistent with what we find in for 
24 < < 171 in the larger masks. From the marginal¬ 

ized distributions we find As < 2 /iK^ at 95% confidence 
and Ad = 0.033 ± 0.008 /iK^, in antenna temperature 
units at the given frequency centers. As it is based on: 
1) Ps and psd from different regions of sky, 2) synchrotron 
and dust power from wide angular scales, and 3) different 
spatial weights than that used by the BICEP2 and Keck 



al Ps Psd 


Fig. 7. — The parameter probability distributions for 50 < ^ < 
no in the of 0.28 region sampled with the MCMC methods. 
The panels on the diagonal show the marginalized distribution for 
each parameter. The parameters at the maxima of the probability 
distributions are shown in blue lines. The off diagonal panels show 
the joint, marginalized distributions, with the contours at 0.5, I.O, 
1.5, and 2.0cr. Note /3s is cut off by the prior distribution. Figure 
made with Foreman-Mackey et al. (2014). 

experiments, it is more of an expectation rather than a 
prediction. Note the dust level is ^2 times smaller in am¬ 
plitude than our = 0.28 region, but it is unclear how 
the synchrotron level scales due to the noise in WMAP. 

6 . DISCUSSION AND CONCLUSIONS 

We have shown a basic polarized foreground model 
consisting of synchrotron, dust, and their correlation 
works for all possible single- and cross-frequency spec¬ 
tra formed with the WMAP9 and Planck 353 GHz data. 
The model needs four parameters: 1) synchrotron ampli¬ 
tude, 2) dust amplitude, 3) synchrotron spectral index, 
4) synchrotron-dust correlation coefficient as a free pa¬ 
rameter with the dust spectrum fixed with the Planck 
dust model, or the dust spectral index as a free parame¬ 
ter with the synchrotron-dust correlation coefficient fixed 
from the definition (equation (7)) using synchrotron and 
dust templates. 

We can draw a number of conclusions from the model. 
1) The foreground minimum for BB in 50 < ^ < 110 for 
the three regions we consider is ^75 GHz, consistent with 
other studies Bennett et al. (2013); Planck Collaboration 
et al. (2015b). 2) There are three primary atmospheric 
windows for observing the CMB at 45, 90, and 150 GHz. 
For the cleanest 50% of sky in polarized dust emission, 
the ratios of the synchrotron to dust amplitudes at the 
three frequencies are approximately 9.0, 0.3, and 0.02 
respectively. 3) In our model, any change in slope of the 
synchrotron index in frequency is explained as a cross 
correlation with dust. Thus, in each Gbin, one index 
describes the synchrotron, although it possibly flattens 
as one moves to smaller angular scales as shown in the 
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Fig. 8 .— Top panels show the best fit foreground model and the uncertainty from MCMC for BB in 50 < £ < 110 for the three masks, 
similar to the bottom right panels in Figures 3, 4, and 5. The model lines follow the same convention as before. Note a difference of 0.5cr 
in (3s for the of 0.28 region can lead to a factor of ~ 4 times larger level of synchrotron power at 100 GHz. The bottom panels show 
the likelihood of the foreground spectrum index between 90 and 150 GHz, along with that of the GMB spectrum at —0.7 in blue. See 
Section 5 for details. 


top panel of Figure 6. 4) The EE/BB ratio of power 
is approximately two although there are large variations 
depending on the sky cut and £ range. 5) When averaged 
over the cleanest 30% of the sky in polarized dust, the 
total foreground emission power at the minimum in 50 < 
£ < 110 is « 0.01 /iK^, corresponding to r « 0.15. 

While the model is an excellent fit to the data, we 
note that foreground emission is subtle. At this point, 
there is no need for an additional component but that is 
simply due to the large uncertainties. One expects the 
synchrotron and dust spectral indices to vary spatially 
and there also may be magnetic dust emission (Draine 
& Lazarian 1999). The current data on polarized syn¬ 
chrotron emission are not constraining enough for an ac¬ 


curate determination of the precise levels in the clean¬ 
est region of sky. Instruments are being built for better 
measurements of the synchrotron emission (e.g., CLASS 
(Essinger-Hileman et al. 2014)). 
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